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Abstract 

An /ip-adaptive Discontinuous Galerkin Method for electromagnetic wave prop- 
agation phenomena in the time-domain is proposed. It differs from many other 
adaptive algorithms in two fundamental aspects: it aims at reducing the true ap- 
proximation error, i.e., it is not based on residuals or heuristic measures such as 
steep gradients, and it does not involve any tuning parameters. We allow for arbi- 
trary anisotropic refinements in the approximation order p and the mesh step size h 
regardless of the resulting level of hanging nodes and, hence, eliminate the neces- 
sity of performing constrained refinements for restoring mesh regularity. This leads 
to meshes with a minimal number of degrees of freedom. The adaptation process 
is guided by so-called reference solutions 1 1 2 1, which are employed for estimat- 
ing the solution error and finding the most suitable type of refinement. During 
mesh adaptation the numerical solution is transferred to the new discretization by 
means of orthogonal projections between Finite Element Spaces. The projections 
preserve the numerical stability of the scheme. Numerical examples are presented 
showing that the algorithm is able of respecting an a priori user-set error tolerance 
throughout time-domain simulations. 

Keywords 

Discontinuous Galerkin Method, dynamical /ip-adaptivity, error estimation, time-domain 
electromagnetics 

1 Introduction 

In this article, we are concerned with solving the Maxwell equations for electromag- 
netic fields with arbitrary time dependence in a three-dimensional domain C M'^. In 
order to achieve this goal the Discontinuous Galerkin Method (DGM) |[3||4j is applied 

*S.M. Schnepp is with the Laboratory for Electromagnetic Fields and Microwave Electronics, ETH 
Ziirich, Gloriastrasse 35, 8092 Zurich, Switzerland e-mail: schnepps@ethz.ch. Manuscript submitted to 
IEEE TAP on Jan 29, 2013. 



1 



on hp-refmed meshes, which dynamically and autonomously adapt as the electromag- 
netic fields evolve. Here, h denotes the (local) mesh step size and p the (local) approx- 
imation order The mesh refinement is driven by a robust local error estimate based on 
so-called reference solutions 1 1 , 2 1 . 

The DG method is a Finite Element type method, which has gained wide acceptance 
as a high order numerical method, which is very well-suited for time-domain problems. 
It combines the usually opposing key features of high order accuracy and flexibility. Its 
flexibility stems from the fully localized character of the numerical approximation. In 
particular, the method can easily deal with meshes containing hanging nodes as stated 
in |]5j, which makes it particularly well suited for ft,p-adaptivity. 

There is a well established body of literature on the DG method for various types 
of problems available. It has been thoroughly investigated by several research groups 
(see e.g. ||5}|7| and references therein). Concerning Maxwell's equations in time- 
domain, the DGM has been studied in particular in [7^10J. The latter two make use 
of hexahedral meshes, which allow for a computationally more efficient implementa- 
tion ITT). 

This paper focuses on error controlled dynamic /ip-adaptation. Mesh refinement 
and specifically /ip-adaptation has received considerable and continuous attention in 
applied mathematics. Despite that, there has not been a lot of development in the 
engineering community. Though a number of works on adaptive meshes with the Fi- 
nite Volume Method are reported in the field of mechanical engineering, especially in 
fluid and gas dynamics (see e.g. p2]-[l4|) publications on /ip-adaptive Finite Element 
Methods with focus on engineering are rare. Published works typically originate in 
the mathematics community. Many of them have a clear focus on the mathematically 
rigorous derivation of error estimates and computable error bounds. 

The first published work on h-, p- and /ip-adaptivity within the DG framework is 
presumably | fT5| , where the authors considered linear scalar hyperbolic conservation 
laws in two dimensional space. Hyperbolic problems have also been addressed, e.g., 
by Flaherty, Shephard and co-workers. In 1 1^ two-dimensional hexahedral meshes are 
considered, where refinement occurs isotropically in h and p. In a follow-up work [fT7|, 
three-dimensional settings were considered. However, in this case the adaptation was 
limited to /i-refinement. In |18 19 1 the DGM was employed. The former contribu- 
tion deals with problems in two dimensions with limited /ip-adaptation, the latter one 
concentrates on two-dimensional problems as well but includes one example using a 
tetrahedral mesh for an essentially two-dimensional forward facing step problem. Re- 
finement is limited to pure /i-refinement. 

A large number of contributions has been authored by Houston and various co- 
workers. They present a number of approaches to adaptivity and deal with first-order 
hyperbolic problems in |20 21), using adjoint solutions |21 22) or estimating errors in 
an energy norm ]23p4[ . The contributions have a clear focus on the rigorous derivation 
of error estimates and error bounds. Applications are limited to one or two space 
dimensions. 

Recently, Solin and co-workers published papers, where they apply dynamical hp- 
meshes for various coupled problems including electromagnetics in two space dimen- 
sions 1 25 26). They employ the concept of reference solutions for controlling mesh 
adaptivity and perform refinements, which are fully anisotropic in both mesh parame- 
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ters h and p. 

To the best of our knowledge, this is the first work presenting error-driven, dynam- 
ical /ip-adaptive meshes for time-dependent, real-world problems in three-dimensional 
space. Mesh adaptation is performed such that a prescribed error tolerance is respected 
throughout the complete simulation. Apart from large savings in computing time and 
memory, the availability of an error estimate is a significant benefit, as it allows for 
reliable statements on the accuracy of the numerical solution. 

This article builds on our previous work published in |27 1, where a general formula- 
tion of the DGM on non-regular hexahedral meshes was introduced. In particular, tech- 
niques for the efficient computation involving non-conforming element interfaces and 
the projection of DG approximations during mesh adaptation were presented. These 
techniques will be applied below. We will not repeat the details but summarize the 
main ideas and steps where we consider it to be helpful. 

We wish to stress that efficiency is not a secondary concern. During a time-domain 
simulation a very large number of element adaptations has to be carried out. Each 
of them is associated with a projection of the local approximation to another finite 
element space (FES). If realistic three-dimensional applications are to be solved, it 
is essential to apply efficient projections techniques, i.e., avoid run-time quadratures. 
The same applies for the computation of surface integrals involving non-conforming 
element interfaces due to hanging nodes. 

In this paper, we propose a robust technique for estimating local and global er- 
rors during time-domain simulations, which is employed for driving the dynamic mesh 
adaptation. The proposed algorithm reduces the true approximation error, i.e., it is not 
based on residuals or heuristic measures such as steep gradients. A second important 
property is that it is entirely devoid of any tuning parameters. The adaptation can be 
performed in four major modes: isotropic in h and p, anisotropic in one of h or p, 
and fully anisotropic in h and p. Unconstrained refinement in h is possible because 
we allow for high level hanging nodes. The number of degrees of freedom (DoF) in a 
discretization will usually decrease from the former to the latter mode, while the com- 
putational load for finding the adapted mesh increases. However, we will show below 
that great savings in both, the number of DoF and computational time can be achieved 
by using fully anisotropic adaptivity. 

The remainder of this article is organized as follows. In Sec. |2] the notation and 
Finite Element Spaces (FES) are introduced, which are applied for obtaining a weak 
DG formulation of Maxwell's equations. Section|3]is devoted to the mesh refinement 
algorithm. First the individual steps, which constitute an adaptive algorithm are dis- 
cussed. They are error estimation, element marking, the /i-p-decision and the actual 
mesh adaptation. For each step a brief description with a review of existing techniques 
is provided, before we proceed with the details of our realization of each step in the 
Sections [3. 2| to [33] Examples are presented in Sec.|4] which include a waveguide and 
an antenna problem. Section |5] summarizes the findings and concludes the article. 
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2 Discretization of Maxwell's Equations 



In the following we assume resting, heterogeneous, linear, isotropic, non-dispersive 
and time-independent materials. Then, the magnetic permeability, /i, and dielectric 
permittivity, e, are scalar values depending on the spatial position only. Under these 
assumptions Maxwell's equations read 

VxE(x,t) = -Ai(x)^H(x,i), (1) 

VxH(x,i) = e(x)^E(x,0+J(x,t), (2) 

with the spatial variable x e and the temporal variable t G [to,T] C M subject 
to boundary conditions specified at the domain boundary dfl and irutial conditions 
specified at time to- The electric and magnetic field vectors are denoted by E and H, J 
denotes the electric current density. We also introduce the electromagnetic energy W 
contained in a volume V given by 



W{t) ^ J^^ (e(x)E(x, tf + m(x)H(x, tf) d^x. 



(3) 



Discretizations of Maxwell's equations using the Discontinuous Galerkin Method 
have been obtained among others in |7j|9][T0]|28|, where the former two employ tetra- 
hedral meshes and interpolatory basis functions. The latter two employ hexahedral 
meshes and interpolatory or modal basis functions, respectively. We will follow the 



framework and notation described in our previous work |27 1, which makes use of hex 



ahedral meshes and modal basis functions as introduced in IITO^ 



2.1 Notation 

We denote by 7/i a tessellation of the domain of interest composed from non- 
overlapping hexahedra 71 such that Tii = UiLi ^ covers Q. The tessellation is re- 
quired to be derivable from a regular root tessellation To by means of elements bisec- 
tions. However, we do not demand the resulting tessellation to be regular, i.e., we allow 
for hanging nodes and specifically for high level hanging nodes. The number of bisec- 
tions performed for obtaining element % is denoted by Li in the isotropic and L^^i in 
the anisotropic case where d corresponds to any of the spatial coordinates {x, y, z}. 

We call the intersection of two neighboring elements Ti OTk their interface lik- 
As we allow for hanging nodes, every face J^j of a hexahedral element may be parti- 
tioned into several interfaces depending on the number of neighbors K such that J-j = 
Ufc=i-^i'=- ™ important difference to most other works including |7[|9||28|, 

which require one-to-one neighborhood relations. The (inter-)face orientation is de- 
scribed by the outward pointing unitary normal . The union of all faces is denoted 
by J-. The volume and edge length measures of element i are denoted by 1 7i | and 1 7d,i | • 

Allowing for high level hanging nodes is very important in the context of adaptivity 
as it allows for unconstrained refinement. 
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2.2 Finite Element Spaces and Approximations 

An essential characteristic of DG methods is that trial and test functions are defined 
with element-wise compact support 

I 0, otherwise. 
Cartesian grids allow for employing tensor product basis functions of the form 

^P(x)= (g) ^P^id), (5) 

d tE {x.y,z} 

where p is a multi-index obtained from all pd = O-.Pd- The local finite element spaces 
(FES) V^(7i) spanned by the basis functions are given by the tensor product of the 
respective one-dimensional spaces 

( V^)r. = (Vf )r.,. ® iK' K..® (Vf ^ )r.,. , where (6) 
(Vf )r.,. - span{v,f (d); < < F^}. (7) 

The approximation may, thus, make use of different orders in each of the coordinate 
directions, where we drop the subscript if they are equal. We do not choose an inter- 
polatory basis but follow a spectral approach and apply Legendre polynomials scaled 
such that (To) 

Pd / \ Qd I \ A \ \ I d.i\i Pd Q_d 

ni.^mi.^)^x^\^ , (8) 

7-^. 10, otherwise. 

Associating an FES (|6]l with each element Ti of the tesselation defines the Finite 
Element discretization, where the electric and magnetic field approximations E/j and 
H/j are represented as 

N 



E(x,t) « E,,(x,t) =0E,(x,O, (9) 

1=1 

N 

H(x,t) « H„(x,t)=0H,(x,t), (10) 



with the element local representations 

E,(x,t) = 5]ef(t)^f(x), (11) 

p 

H,(x,t) = ^hni)^r(x). (12) 



The time-dependent vectors of coefficients e = (e^, .., ef , .., e^J^, .., e^)^ and h 
(h^, .., hf , .., h^l^, .., h^)^ are the numerical degrees of freedom. 
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2.3 Weak DG formulation 



We are now in the position to state the weak DG formulation of Maxwell's equations. 
Following the Galerkin procedure ([TJ and (|2| are multiplied by a test function tp and 
integrated over the domain Vl. Due to the compact support property Q the integration 
can be carried out over every element Ti individually. Next, we perform integration by 
parts of the curl-terms and replace the exact field solution for the approximations (|9]l 
and ( [TO] i. This leads to the semi-discrete variational problem of finding e and h such 
that 

j^^ i; /i^H^ d^x - J^^ ( W) X d^x 

+ / V(nxE,Od"x = 0, (13) 

■JdTi 



d C 
V e^E,, d^x + j^^ (W^) X d^x 



/ V(nx H^)d2x = 0, (14) 

JdT, 



Vi = 1, .., iV; V?A e Vi. For the above equations to be well-defined it is required that 
ip E in the interior of %, which is fulfilled for the chosen Legendre basis. Note 
that E/i and H/, denote the numerical trace of the electric and magnetic field, which 
is single-valued for each vector field component at element boundaries. Introducing 
the numerical trace is a necessary step for resolving the ambiguity of the numerical 
approximations (|9]l and ( [TO] i at element interfaces. Due to the definition of the basis 
function support in Q, the components for the vector fields E;i and Hh are single 
valued at all points x e T\J- but double-valued for all x e J^. The numerical trace is 
computed as 

^ {{yE}}x., n,fc X |HIx., 

Typical choices are the centered and upwind value obtained by setting 7 to zero or one, 
respectively, where the upwind value is the solution of the Riemannian problem p9) . 
Above {{ •}} and | ] denote the average and jump operators 

{{a}}i., = (afeii^, + a,|i, J/2, (17) 
Ni., = afe|i^^ - a^ix.^. (18) 

The intrinsic impedance and admittance are given as 
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The surface integrals in ( [T3] l and ([T4| represent interelement fluxes, the volume in- 
tegrals are referred to as the mass and stiffness terms according to standard FE nomen- 
clature. In the following the dependence of the spatial and temporal variable is not 
written down explicitly. 

Note that no assumptions on the grid regularity have been made in the derivation. 
This is in a sharp contrast with Finite Element Methods based on edge elements, which 
require augmentation by edge constraints if hanging nodes are to be included ||2]. In 
DG-type methods non-regular grids are no methodological issue, they only make the 
implementation more involved. However, in p7| , we proposed techniques for dealing 
with high level hanging nodes in an efficient manner. The relative ease of handling 
non-regular meshes combined with the strictly element-local character of the numerical 
approximation make DG methods an ideal candidate for /ip-adaptivity. 

3 Automatic and dynamic /ij9-adaptation 
3.1 Steps of the Adaptive Algorithm 

Devising an /ip-adaptive algorithm requires four major steps. 

1 . Derivation of global and local error estimates 

2. Definition of a marking strategy for assigning a refinement/derefinement label to 
each element 

3. Deriving criteria for making the /i-p-decision 

4. Definition of the actual mesh refinement/derefinement operators 

For each of these steps several alternatives are possible. We will briefly list a few 
popular techniques and describe the main underlying idea before naming the approach 
followed in this contribution along with the reasoning behind this choice. 



3.1.1 Error Estimation 

Error estimators or indicators can be obtained by expressing a residual through the 
numerical approximation. Residual based estimators in the context of Maxwell's equa- 
tions have been developed, e.g., in |23||30J for Nedelec Finite Elements and DG re- 



spective, in 1 3 1 32 1 in the context of DG methods with applications outside electrody- 



namics, or recently in |33| and 1 27 1. In the latter contribution, we adapted the residual 



estimator of 131] to Maxwell's equations. 

Highly accurate estimators can be constructed based on adjoint solutions |34j35). 



where the latter one is applied in a DG setting. However, the accuracy of adjoint based 
estimators comes at the price of having to repeatedly solve for the adjoint problem in 
addition. 

This list is far from being exhaustive. More comprehensive overviews are found 
in IIIHID. 

In this article we employ the concept of reference solutions |[T][2][39) for obtaining 
error estimates. A reference solution is a numerically computed approximation, which 
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is assumed to be significantly more accurate than the present approximation. This 
can be achieved by performing one isotropic /i-refinement combined with increasing 
the approximation order by one on the element under consideration. We apply the 
concept in its original form for finding an initial hp-mesh and propose a modified, 
computationally much cheaper variant, which is applied during the transient analysis. 
We found this estimator to be very robust and find reliable estimates independent of 
the local solution smoothness. This is an important advantage over the residual based 
estimate proposed in |27J . 



3.1.2 Marking Strategy 

Once the element-wise approximation error is estimated, each element is assigned one 
of the labels refine, derefine or retain according to the marking strategy. The marking 
strategy, hence, has a strong impact on the number of DoF in the computational mesh. 

Popular strategies include error equidistribution, the fixed fraction strategy or vari- 
able fraction strategies such as bulk-chasing, commonly known as Dorfler-marking. 

The goal of the former strategy is to equilibrate the local errors by refining or dere- 
fining elements such that « TOL/V^iV, where Si is the local error estimate and 
TOL is a user-defined error tolerance f40l. For the fixed and variable fraction strate- 
gies, the elements are ordered by their estimated error at each refinement step. Then, 
for the former approach, a fixed fraction of elements form the top and bottom are 
marked for refinement and derefinement. The variable fraction or Dorfler-marking on 
the other hand continues to mark elements from the top and bottom of the list until 
their accumulated error accounts for a certain percentage of the total error This can 
be expressed as finding a minimal subset 7^+ and a maximal subset Tf^ of Th such 
that X]7-g7-{+'-> ^1 — ^{+ -} I^TiGTh ^i' ^here the sign indicates refinement and 
derefinement. As the values of 0{+.-} indicate fractions of the total error, the Dorfler- 
marking can be considered as a fixed fraction marking with respect to the total error. 

Often a few percent of the elements make up for more than 90 % of the total error, 
while most of the elements contribute to the total error by less than 5 %. As the situa- 
tion also might change throughout a time-domain simulation, we consider the variable 
fraction marking the most suitable for our problems. 



3.1.3 The /ip-Decision 

Following the decision on which elements to adapt, the kind of adaptation has to be 
chosen, i.e., h- or p-adaptation. This decision is guided by the local solution smooth- 
ness. It is well known that for sufficiently smooth solutions consecutive p-enrichment 
leads to exponential convergence, whereas ^.-refinement yields algebraic convergence 
rates only |41 42]. This is reflected by an a priori error estimate given in presumably 



the first paper on /ip-adaptivity in the DG context [ 15J and others. We repeat it here in 
the simplified form 



MUp < C i^Jjrj II^IU,,;, (20) 

where C is a constant independent of h and p, ^li = mm{pi + l,s),i/ = s — 1 with s the 
Sobolev regularity index, u is the exact solution, and ||| • ||| indicates a mesh dependent 
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Figure 1: Global projection error of a Gaussian and a trapezoidal waveform as depicted 
in the insets. The plot in the top panel shows convergence of the error at a rate of P + 1 
for the Gaussian waveform. In the bottom panel, the low regularity limits convergence 
to first order. 

norm. In this norm the mesh parameters h and p vanish if the continuous case is ap- 
proached, i.e., in the limit of /i — 0, p — ?► oo. The above a priori estimate establishes 
convergence rates of the method and clearly predicts how the error in the numerical 
solution behaves under h- or p-refinement depending on the local solution regularity. 
Specifically, it states that the achievable convergence rate is limited by the local reg- 
ularity. Raising p does not yield better convergence if the solution is not sufficiently 
regular. In fact, it often leads to problems in the form of oscillations of the numerical 
solution. 

FigurefTlillustrates the dependence of the convergence rate on the regularity as pre- 
dicted by ( |20l ). The waveforms depicted in the insets, i.e., a Gaussian and a trapezoidal 
waveform in one-dimensional space, are projected to spaces with P varying from 
zero to five. The plots show the global error measured in the L^-norm. While the 
convergence rate increases from one to six with every increase of P for the Gaussian 
waveform, convergence is limited to first order in the latter case. Both results agree 
with the prediction obtained from ( pOj i. We note that the estimate ( pOj i cannot be ap- 
plied for assessing the accuracy of a given numerical solution since it involves unknown 
constants and the exact solution. 

However, exponential convergence in terms of DoF can be obtained even for lo- 
cally non-smooth solutions as well by employing proper /ip-refinement |41 1. To this 
end, regions of low regularity are embedded into /i -refined areas of the mesh using low 
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order polynomials. This encapsulates the non-smooth regions and limits their effect 
on the global convergence. Then, p-refinement is applied everywhere else. Thus, the 
performance of the adaptive method critically hinges on correct /ip-decisions. In or- 
der to be in the position of performing anisotropic /ip-refinement in three-dimensional 
space, we require information about the directional smoothness of the unknown solu- 
tion. This means that unlike for error estimation, more than one value per element is 
required. This is the main reason, which renders the task of smoothness estimation 
more difficult than error estimation. 

A variety of techniques have been proposed for assessing the local smoothness or, 
more general, for making the /ip-decision. The simplest ones makes use of information 
that is available a priori such as the position of field singularities due to corners and 
spikes | |43| . However, instead of relying on geometric information, we rather wish to 
drive the /ip-decision based on the actual numerical solution. 

Known methods include the type parameter technique | |44J , 'Texas 3-Step' l|45j, 
mesh optimization techniques | [46) , error prediction pT] or local regularity estima- 
tions |48 -51 1. Descriptions of all methods are beyond the scope of this paper, and we 
refer to |52| for an extended overview including descriptions. 

However, a popular method is the estimation of the Sobolev regularity index s in 
a local manner One viable technique for achieving this goal is described briefly in 
the following, as it is illustrative for understanding why we pursue a different strategy. 
We focus on |49 50], where the authors first develop such a strategy based on moni- 
toring the decay rate of the sequence of coefficients in the Legendre series expansion 
of the numerical solution. In the second contribution this strategy is amended by the 
estimation of the analyticity in a certain neighborhood by estimating the local size of 
the Bernstein ellipse. If the solution is found to be analytic in the neighborhood p- 
refinement can be applied and the estimation of the regularity index is skipped. This is 
shown to produce meshes with slightly less DoF for a given error tolerance. The draw- 
back of this method and many others is that a certain number of coefficients is required 
for the computation of the coefficient decay rate to be robust. Taking into account that 
the Legendre coefficient of order zero provides information about the average in the 
element only, coefficients providing actual decay information start with the order of 
one. Hence, decay rate estimations require second order approximations as a minimum 
(as used in | [49) ) although higher order approximations will make the method more 
robust. Problems also occur if the solution exhibits a pronounced odd-even character- 
istic | [53| leading to an alternation of small and large valued coefficients. The extension 
to problems in two- and three-dimensional space is certainly possible but not unique, 
and the technique loses part of its clarity. Moreover, approximation orders of at least 
two have to be applied in all directions. This leads to a significant number of DoF also 
in elements, which do not require it. 

As we wish to employ approximation orders as low as possible (ideally order zero) 
everywhere the solution permits, we follow a different approach. In |27|, we proposed 
a smoothness indicator, which works down to approximation orders of zero. It was 
originally proposed in |54| for problems in fluid dynamics and adapted to Maxwell's 
equations. It is based on the exploitation of a superconvergence property of the DG 
mefliod (55). 

In this contribution, we employ reference solutions for finding the most suitable 



10 



refinement from a list of candidates. This approach circumvents the issue of regularity 
estimation and the associated difficulties by testing various h-, p- and /ip-candidates 
with respect to the reference solution. With this strategy, the best candidate naturally 
arises as the one offering the best ratio of approximation error £c to the logarithm of 
its number of DoF (e^/ log(#DoF)). The logarithm is required as it is the aim to 
achieve exponential convergence. For finding the initial mesh, we follow the original 
algorithm but propose a computationally cheaper version to be used during the time- 
domain simulation. 



3.1.4 Mesh Adaptation 

More traditional approaches to mesh refinement are longest edge bisection, or red- 
green refinement p6| , however, at this point sufficient information is available for per- 
forming error driven /ip-adaptation, which can be anisotropic in both mesh parameters 
h and p. The high degree of localization of the DG method turns mesh adaptation into 
a purely element-local operation. 



3.2 Error Estimator 



As mentioned above, we employ the concept of reference solutions (cf. |T||2 46 1) for 
obtaining element-wise error estimates. The underlying idea is to use a reference mesh 
%ef, which can be obtained by performing one global uniform refinement step in h and 
p. Obtaining a reference mesh by pure p-enrichment has been proposed as well. Both 
techniques provide a reference mesh based on hierarchic FES enrichment. 
The aim then is to find the minimal hp-mesh such that 

\\4n = l|u - u^pIIt^, = ( E II" - u'^pllr. ) ' < TDL, (21) 

where || • Hj-^ denotes the global L^-norm, and it is taken into account that the solution 
is a vector field. In the following u is used for denoting the electromagnetic solution 
(E, H). As only the approximation u^p is known but not the exact solution u the target 
( |2l] l cannot be achieved directly. However, it can be achieved asymptotically as 

\Mn = ( E ll"'-rf - n,pu,ef||^^) ' = ( E ll^'^pllr.) ' ' (22) 

where Hhp is a projection operator from the enriched reference FES Vret to a space Vc 



associated with a refinement candidate, which will be defined in Sec. 3.5 The space 
Vc is reduced with respect to the reference space but enriched with respect to V-Ti such 
that 

Vt. C Vc C Hef. (23) 

The refinement candidates are obtained by refining element Ti in h and/or p such that 
( |23] l holds. Refinement can be anisotropic in one or both of the mesh parameters. 
From ( |22l ) it follows that the element-wise error estimate is given as 

Ik/ipllr. = ||Uref - n^pUrefllTi- (24) 
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Given a reference solution, the global and local error estimates ( p2] l and (|24]| are fully 
computable. 

We pursue different strategies for obtaining a reference solution during the con- 
struction of the initial mesh and the iterative time stepping. This is owed to the fact that 
the initial data, which is required for constructing the initial hp-mesh is known while 
the transient solution obviously is not. 

3.2.1 Initial Mesh 

Starting from the root tesselation To with some uniform polynomial order Pq a refer- 
ence mesh is constructed by performing one isotropic refinement step in h and p. We 
note that this has not to be done globally, but it can rather be done consecutively with 
each element of the current tesselation. The DG approximation /t- of a given function 
/ on the element %, is obtained by applying the orthogonal projection operator 11 

where (u, v)j'. denotes the inner product J!^ uwdx on the element %■ Equipping the 
FES (|6| with an inner product defines a Hilbert space. Hence, the above projector 
yields the best approximation in the L^-sense. After projecting the initial data to the 
refined elements the approximation error £i of element % is estimated using (|24]i. This 
procedure is repeated for all elements of the current tesselation. The global error £hp 
is obtained from ( |22] i. The construction of the initial /ip-mesh terminates when the 
stopping criterion e^p < TQL is met. 

3.2.2 Dynamical Mesh 

In the construction of an optimal initial hp-mesh the reference solution at each iteration 
can be generated because the initial data is known exactly. Obviously, this approach 
cannot be transferred immediately to the transient analysis. In |25l, the authors ap- 
proach the transient case by employing Rothe's method. In contrast to the widely used 
Method of Lines, Rothe's method discretizes the time variable first while preserving 
continuity of the spatial variable. Then, in every time step, the evolutionary PDE is 
approximated by means of one or more time-independent ones. The number of time- 
independent equations per time step is proportional to the order of accuracy of the time 
discretization method. This approach allows for applying the same techniques in the 
transient analysis that were used for obtaining the initial mesh at the cost of having to 
solve for systems of equations in every time step. 

For performance reasons we prefer to employ explicit time-integration. The straight- 
forward extension in this case is to compute on two meshes, the hp-mesh fulfilling the 
error tolerance and its reference mesh. However, this approach is prohibitively expen- 
sive, both in computing time and memory consumption as the reference mesh usually 
has about ten to 60 times more DoF depending on the approximation order Taking 
into account that, moreover, the reference solution largely exceeds the required accu- 
racy and is employed for driving the adaptivity only, we seek a different approach. 
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To this end, we switch roles of the reference mesh and the mesh used for estimating 
the local error To this end, we claim that the approximation on the current hp-mesh 
is sufficiently accurate for serving as the reference solution and estimate the element 
error by comparing to a reduced FES. This FES can be obtained by dere fining the mesh 
in h and p, or in a computationally significantly more efficient manner by reducing the 
approximation order P. The element-wise error estimate is computed as 



k/ipllr, — l|Ui-ef — HpUi-efllT-i- 



(26) 



Here, the solution on the current hp-mesh is the reference solution and Hp is the pro- 
jection operator to the p-reduced FES. Computing the estimate ( |26l ) is very cheap as it 
comes down to considering the highest order terms of the current approximation only 



p p-i 

p=0 J9=0 



(27) 



where denotes the vector of coefficients of order p local to element %■ Recalling that 
p and P are multi-indices as defined in (|5]l, the local error estimate ( |27] i is computed as 
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(28) 



Evaluating the L^-norm and inserting the scaling property (|8]l of the basis functions 
yields the following form of the estimate 
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(29) 



It is an important feature that the computation of the estimate as given in ( [29| ) is highly 
efficient as no runtime quadratures have to be performed. The global estimate reads 



|e/ipllrh 




1/2 



£?ipllri 



(30) 
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We admit that the approach of projecting the solution to a reduced FES instead 
of an enriched one will negatively affect the accuracy of the error estimation. How- 
ever, it drastically reduces computational costs rendering the method applicable for a 
much larger class of real world problems. In Sec. [4] we demonstrate the robustness and 
reliability of this approach. 

3.3 Marking Strategy 

Also for the element marking distinct strategies are applied for constructing the initial 
hp-mesh and during the transient analysis. 

For generating the initial hp-mesh, we perform variable fraction or Dorfler-marking. 
The number of mesh adaptation iterations required for obtaining the initial mesh de- 
pends on the fraction of the total error. Less iterations are performed for large fractions. 
However, this usually leads to a slightly larger number of DoF. For the construction of 
the initial hp-mesh elements are marked for refinement only. 

During the transient analysis a slightly altered marking strategy is employed. This 
strategy is a variable fraction strategy with respect to the number of elements as well 
as to the total error. For every mesh adaptation a minimal subset 7^^ is assembled such 
that 

J2 £'>min( ^ 02^2 max (4p-T0L2,0)). (31) 

r,er+ r^^n 

Loosely spoken, the above equation states that the size of the minimal subset of ele- 
ments to be refined is not larger than determined by the given fraction 9, but it can be 
smaller if the global error is close to the prescribed tolerance. If the estimated global 
error is smaller than TOL, the set is empty, and no elements are refined in this adaptation 
step. If we were to apply the marking strategy in the same way we did for obtaining 
the initial mesh, the algorithm would continue refining elements even if the estimated 
error is less than the tolerance. 

As stated above, we assume that the approximation on the current hp-mesh is suffi- 
ciently accurate for serving as the reference solution. This statement should ideally be 
true for every element. Therefore, marking elements for derefinement has to be done 
with care. We recall that mesh adaptation during the transient analysis is a dynamic 
process. Therefore, elements suitable for derefinement, which are not marked as such 
in an adaptation step are again considered for derefinement in the next step. In the 
examples in Sec.|4] we show that the mesh derefinement works well despite the careful 
approach. 

3.4 The /i-p-Decision 

If a refined reference solution is at hand, the task of performing a regularity estimation 
or a similar task becomes obsolete. Given a tesselation 7/i the best refinement option 
is found by testing a list of candidates against the solution on the reference mesh 77ef • 
The most suitable refinement is the one offering the best ratio Ec/ log(#DoF) with the 
candidate error Ec- 
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The size of the Hst of candidates can vary considerably. It depends on the global 
refinement strategy, i.e. isotropic refinement only, fully anisotropic, or anisotropic in 
one of h and p only, but it also depends on the permissible increment and decrement 
in the /i-refinement level AL and AP. In this paper, we restrict both to one. However, 
candidates have to be competitive. This means that increasing the /i-refinement level 
L, or Ld in the anisotropic case, goes along with a reduction of P in order to prevent 
a strong increase of the number of DoF in the element. The approximation order is 
reduced such that the number of DoF of the candidate is as small as possible but larger 
than the one of the current element (#DoFc > #DoFj). If isotropic refinement is 
applied only, this limits the number of refinement candidates to two, one /i-candidate, 
consisting of eight elements with possibly decreased approximation order P, and one 
p-candidate. For fully anisotropic refinement a number of fourteen candidates is con- 
sidered, which are obtained by refining each of the mesh parameters h or p in one 
direction (three candidates), two directions (three candidates), and all three directions 
(one candidate). The approximation order P of /i-candidates is reduced as described 
above. We do not construct more candidates by allowing for larger AL and AP as 
testing the candidates is computationally rather expensive. 

The procedure above applies to mesh refinement. For the case of mesh derefine- 
ment, it is natural to proceed in a similar manner and set up derefinement candidates 
with a smaller number of DoF. 

In the dynamic case, the procedure requires modification as the problem is en- 
countered that a refined reference solution cannot be constructed. An error estimate is 
obtained by projecting to a p-reduced FES, however, given our description of regularity 
estimation based on coefficient decay rates in Sec. |3.1.3| it is doubtful that regularity 
information can be extracted from a comparison of two solutions of the order P and 
P — 1 in a robust way. As there is not enough information available for making a rea- 
soned /ip-decision, the respective element is refined uniformly in h and p in order to be 
on the safe side. During the next mesh adaptation one or more of the refined elements 
might be derefined again according to the estimated error. Derefinement can be carried 
out anisotropically depending on the global refinement strategy. During the transient 
analysis anisotropic adaptation, hence, occurs during derefinement only. 

3.5 Refinement and Derefinement Operators 

Upon mesh adaptation the numerical approximation given on the current ftp-mesh Th 
has to be transferred to the adapted mesh Th*- The objective is to find the best rep- 
resentation of Ufip on Th* with respect to the i^-norm. For all adaptations (h/p re- 
fine/derefine) this is achieved by applying the orthogonal projection operator 11 intro- 
duced in p5| ). Due to the compact support of the basis, an unconstrained projection 
can be carried out in a strictly element-wise fashion. Additionally, the tensor product 
property of the basis (|6]l allows for performing the projection along each dimension 
individually. This reduces the three-dimensional quadrature of complexity order three 
in the number of quadrature nodes to a product of three one-dimensional quadratures 
of complexity order one. 
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3.5.1 /i-Adaptation 

The one-dimensional FES associated with the left hand and right hand refined element 
Td,i,Td,r are denoted as Vd,i, Vd,r with V+ — Vd,i ® Vd,r being their direct sum. The 
approximations on the refined elements of % are obtained as 

{uHpYr, = (tfu„p)^, (u„p)^ = (Tl'nhph, (32) 

with the orders I < Ld,r < Rd and Ld, Rd the maximum approximation orders in 71 
and %■ Due to the tensor product character of the basis, this can be expressed as 



7T 



p 

for the left and right child, respectively. The projection operator reads 

Using the orthogonality of the basis ([H), the projection in one coordinate, e.g. x, sim- 
plifies to 



Px 



Kpv Sl^P^ H < ('^l" ' "^^i Tx , ' (3^) 

' Px 



where (5 denotes the Kronecker delta. Note that above the summation index is px 



whereas in (33 i it is p. This reduces the number of addends from \P\ to P^. and sig- 
nificantly reduces the computational costs. As {ip\'' , ip^'' , is identically zero for any 
Px < Ix, the above sum can be limited further to the range [Ix, Px], which reduces the 
number of addends to the minimum possible. 

For the merging of elements, the approximation within the parent element, %, is 
given as a piece-wise defined function within its child elements. Despite that the func- 
tion to be projected is not continuous, the projection operator 11 is applied in a fully 
similar manner. The projection reads 

(u,p)^^ = (nfu,p)r. = (nf (u„p)^)^_ + (n^'(u;,p)^)^_ 

^ ^u((nv,')7i + E<(nvn7; (36) 



where the simplifications (33 i and (35 i apply. 

If the approximation orders L and R of the refined elements are such that L, R> P 
the FES (V+)ri is a superspace of (V^)^;- Then, the projection is point-wise ex- 
act in the /i-refinement case. Point-wise exactness is in general not obtained for h- 
derefinement. Apart from the FES nesting argument, this is also clear from the fact 
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that the approximation is allowed to be discontinuous across the interface of the re- 
fined elements, but it has to be continuous after merging. 



3.5.2 p-Adaptation 

The adaptation of the local approximation order P is achieved by enriching or reducing 
the local FES followed by a projection of the current approximation to the adapted FES. 
Enriched and reduced spaces are obtained as 

Vf+i = V,ru{^r^}, (37) 

vf-' - vr\wr}. (38) 

The tensor structure also allows for enrichments and reductions in one or two spatial 
dimensions only. 

The basis Q-l^]) is constructed from Legendre polynomials of increasing polyno- 
mial orders. The orthogonality property (|8]l implies that the basis is also hierarchic, 



which largely simplifies p-adaptation. Formally, we perform the projection ( 25 1 but 
by virtue of the orthogonality, this equals extending the vector of local coefficients 
Ui with the new coefficients u^"*"^ during enrichment and removing the highest order 
coefficients upon reduction. The new coefficients are initialized to zero. All existing 
coefficients require no modification, as they are unaltered under projection from the 
current to the adapted FES . 



3.5.3 Comments on Practical Issues 

During one time-domain simulation a very large number of adaptations is performed. 
These have to be administered in a way, which allows for an efficient traversing of 
all elements in each time step. Additionally, parent-child information is required for 
simplifying mesh derefinement. In this context, tree structures emerge as a suitable 
storage format. They allow for operating on the current discretization by working on 
the tree leaves only but contain the refinement history and parental relationships as 
well. 

In the case of isotropic ft -refinement, the tree is organized using an octree-structure, 
where each of the eight children is assigned to one branch. In an octree-structure, every 
element and its associated node is either a non-reducible element of the root tesselation 
To or one of eight children of a single parent element. The depth of a node in the tree, 
i.e. the number of ancestor elements to the respective root element, corresponds to the 
number of consecutive /i-refinements. This has been defined as the ^.-refinement level 
L before. 

This organized view breaks down if non-anisotropic refinement is permitted. We 
refer to Fig. |2] for the following explanation. For the sake of clarity, a single two- 
dimensional element is considered. In Option I of the left hand side example only 
anisotropic /i -refinement is applied. We extend the mesh representation tree in the same 
way as for isotropic refinement, i.e., the splitting of elements for every /i-refinement 
is represented by extending the tree downwards from the respective node. In Option 
II, a combination of anisotropic and isotropic ft-refinement is performed. This yields 
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the same final mesh but a different representation tree. For the example on the right of 
Fig.[2j the same refinements are performed in a different order leading to identical final 
results and apparently identical representation trees. 

The problem associated with the representation trees in Fig. [2] becomes visible 
when we attempt to derefine the mesh. This is achieved by cutting branches from 
the leaves upwards to the root, which immediately implies that derefinement has to 
occur exactly in the reversed refinement order. This is not a desired behavior, as the 
solution can develop in a way such that a different derefinement order would be more 
suitable. We wish to point out, that the simple two-dimensional examples of Fig. [2] 
suggest that this is a minor issue. Nevertheless, in more complex situations in three- 
dimensional space it is a clear disadvantage if mesh refinement and derefinement have 
to be performed in reversed order. 

The issue can be faced in a number of ways, where many of them are computa- 
tionally expensive. As one example, graph theory could be applied for generating a 
new minimal representation tree after each refinement. Our approach is computation- 
ally much cheaper and aims at constructing representation trees of minimal depth. The 
idea is illustrated in Fig. [3] In order to obtain a tree of minimal depth, a new gener- 
ation of children is spawned only if the maximum /i-refinement level L — max(Lrf) 
is increased. For the minimal tree I (Min. Tree I), this is the case for the first two 
refinements but not for the last refinement step. Initially, elements are twice refined 
horizontally yielding a maximum ^.-refinement level of L = 2. The last refinement oc- 
curs vertically. As this does not cause L to increase, the respective leaves are connected 
with the same parent node. This strategy yields identical minimal representation trees 
I and II. However, the uniqueness of minimal trees is not guaranteed by the approach 
as demonstrated in the right hand example of Fig. [3] Nevertheless, for general refine- 
ments in three dimensions trees of a significantly smaller depth are obtained. They also 
provide a more intuitive representation as the tree depth connects with the maximum 
/i-refinement level. 

The true benefit of constructing minimal trees becomes visible when mesh dere- 
finement is considered. In contrast to the trees constructed in Fig. |2] the derefinement 
order is not strictly prescribed by the refinement order The minimal tree I allows for 
derefining such that the meshes at steps one or two are obtained. Additionally, a mesh 
with one horizontal and one vertical refinement of the right hand side element is ob- 
tained naturally. Using minimal trees, identical representations, such as I and II, always 
offer identical derefinement options, which is a significant advantage regarding the im- 
plementation in a computer code. Considering minimal tree III all meshes depicted in 
either of option I or II can be obtained by derefinement. Additionally, other meshes are 
possible by performing derefinement of some elements only. 

The selection algorithm for the most suitable derefinement candidate is depicted 
in Fig. |4] where the mesh of the right hand example in Fig. [3] and minimal tree III is 
considered. Only the right hand half is depicted as derefinement of the other half is 
carried out analogously. Given the current discretization and its representation tree, we 
move one level upwards in order to obtain the topological parent element. The parent 
is the first /i-derefinement candidate. Then, successively all possible /i-refinements of 
the parent are performed such that Ld.c < Ld is respected. The additional ^.-candidates 
for the considered example are depicted in the third row of Fig. |4] Figure |4] depicts an 
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example only. The number of /i-candidates is not fixed and depends on the number of 
/i-refinements of the current discretization with respect to its topological parent. In a 
third step, the purely topological /i-candidates are assigned FES of different orders P 
yielding /ip-candidates. We restrict the generation of /ip-candidates in the sense that all 
elements of a candidate have the same order P. However, each /i-candidate has its own 
P dictated by the requirement that the number of DoF of the candidate is smaller than 
in the current mesh. In a last step, we compute Ec/ log(7^DoF) for each /ip-candidate 
and choose the best derefinement option. 

We note that the above assumes that all elements have been marked for derefine- 
ment. Currently, /i-derefinement is performed only if all children are marked for dere- 
finement. More sophisticated coarsening strategies are subject of current work. 



3.6 Efficient Computations on Adaptive Meshes 

DeaUng with adaptive meshes containing hanging nodes in an efficient way is discussed 
in details in f27 |. In this section, we summarize the most important results only. 

In every adaptation step a large number of projections of the kind ( (32] l and ( [34] l 
has to be carried out. Each of these projections is associated with the evaluation of a 
number of the inner products of the kind 

(^h^r')r.,,' K^'/'^)r./ (39) 

The evaluation of these integrals in a computer code at runtime should be avoided as 
quadratures are computationally very expensive. However, the projection operator for 
performing /i-adaptation is fully determined by the basis given on the current element 
and the adapted elements. As it is independent of the actual solution, it can be ana- 
lytically precomputed for all combinations of basis functions yielding matrix operators 
n. This eliminates the necessity of runtime quadratures. Projections to adapted ele- 
ments reduce to the evaluation of matrix- vector products of type IlUi. This procedure 
respects the electromagnetic energy (|3]l of the current field solution as a strict upper 
limit and was shown to preserve stability of the employed DG time-domain method. 

In a similar way, the interface integrals occurring at non-conformingly refined 
faces, i.e. faces containing hanging nodes, can be precomputed. This levels the costs 
of the flux evaluation at conforming and non-conforming interfaces. 



In Sec. 3.5.2 we pointed out that p-adaptation is a trivial task from the mathe- 
matical point of view. However, its efficient implementation in a computer code is an 
issue. Extending a local vector of coefficients u,; requires storage for the additional 
DoF. Straightforward solutions to this issue are memory resizing through system com- 
mands or the allocation of sufficiently large memory chunks upon element generation. 
The former option is very inefficient as memory resizing through the operating sys- 
tem is slow, the latter one is prohibitively expensive regarding memory consumption. 
We resorted to handling all dynamic memory allocations in our implementation p7] 
through a library based on memory blocking | |58| . See p9] for more details on the 
latter subject. 
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4 Examples 



4.1 Fundamental mode in a waveguide 

As a first example we consider the propagation of a wave packet in a rectangular 
waveguide. We consider a waveguide of type WR 19 working in U-Band. The cutoff 
frequency of the fundamental mode is 40 GHz. The frequency limit for single-mode 
operation is 60 GHz, and the wave packet considered has a frequency range of 45-59 
GHz. The waveguide aperture dimensions are 4.78 x 2.39 mm, and we consider a total 
length of 1 m corresponding to approximately 170 wavelengths. Through this rather 
academic example we wish to demonstrate that the proposed algorithm can cope with 
situations where a very large number of adaptations has to be performed. Throughout 
the simulation the error tolerance has to be respected. Also, the number of DoF should 
remain approximately constant as the wave packet will largely keep its shape. 

The generation of the initial hp-mesh requires 28 iterations with the fraction 9 as 
described in ( (3T] i set to 0.5. The series depicted in Fig. |5] shows the hp-mesh and the 
respective approximation of the -component on the uniform root tesselation, at an 
intermediate iteration and the final hp-mesh. Refinement is allowed to be anisotropic 
in both mesh parameters, h and p, though the algorithm applies no /i-refinement in this 
case. This is reasonable as the solution is smooth. For depicting anisotropic hp-meshes 
we make use of a common visualization technique |jT]j2]. To this end, each face is split 
into four triangles. The tensor product orders are coded with the triangle color If the 
base edge of a green triangle is aligned with the x-axis, then = 4 according to the 
color legend. In the same way = 5 is represented with an orange triangle having 
its base edge aligned with the z-axis. This visualization allows for representing the 
orders in one plot and also gives an immediate impression of predominant directions 
regarding the approximation orders. 

The highest orders in the initial mesh are Px ^ 5 and Pz = 6. As the fundamental 
mode shows no variation in y-direction, no increase of Py occurs. The construction of 
the initial mesh requires seven seconds and yields close to 135,000 DoF. If we allow 
for isotropic refinement only, an initial mesh with 285,600 DoF is obtained within 
twelve seconds. Figure [6] shows the convergence graph of the approximation error 
with the number of DoF in a semi-logarithmic plot, where each circle represents one 
iteration. In this graph, the error reduction occurs along an almost straight line showing 
exponential convergence, which is the desired result. 

Next, the time-domain simulation is performed. Figure|7]shows the _Ey-component 
and the hp-mesh after the packet has traveled to the center and to the end of the waveg- 
uide. The performance of the adaptive algorithm is illustrated more quantitatively in 
Fig.js] The top plot shows the evolution of the estimated global L^-error normalized to 
the error obtained on the initial mesh. The middle and bottom plot depict the number 
of elements and DoF throughout the simulation. The data corresponds to 50 samples 
in time. The dispersion, which can be observed in Fig. |7]is a physical effect due to 
waveguide dispersion, not a numerical artifact. 

In order to assess the reduction in computing time and memory consumption due to 
adaptivity, simulations on a fixed mesh were carried out. The choice of the approxima- 
tion orders was based on the initial hp-mesh depicted at the bottom of Fig.|5] where the 
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Table 1 : Performance of simulations of example 4. 1 using fixed and adaptive meshes 



Orders 


DoF 


Memory 


norm. 


L^-error 






/MB 


Runtime 


/IQ-^ 


3/3/3 


808 


30.4 


2.7 


15.7 


5/1/6 


1131 


35.5 


4.3 


0.13 


5/5/5 


2911 


62.7 


10.3 


1.36 


61616 


4620 


88.7 


20.6 


0.13 


hp 


125-140 


4.7-5.3 


1 


1.01 



highest orders are (5/1/6). The number of DoF, memory consumption, runtime and 
error estimates after the final time step for various fixed meshes are listed in Tab.[T] Us- 
ing approximation orders (4/1/5) globally yields an estimated error approximately 30 
% larger than the /ip-solution, whereas the error using orders (5/1/6) is considerably 
smaller. Regarding runtime the fixed mesh solutions require about three, respectively 
four times longer than the adaptive solution. Memory consumption is higher by factors 
of about six and eight, respectively. This has to be put in relation with the part of the 
mesh that is refined. The length of the refined region is approximately 0.08 m, i.e. 8 
% of the waveguide length. Code profiling showed that about 15 % of the computing 
time is spent for adaptation related tasks. 

For comparison, we also included two simulations using uniform orders of five and 
six. As Tab.[T]shows this does not decrease the error but leads to a significant increase 
of the number of DoF and runtime. This example considering the fundamental mode 
clearly favors anisotropic orders. Nevertheless, anisotropic orders and anisotropic re- 
finement has a big potential regarding savings in computational resources in most ap- 
plications. This is consistent with the findings of |^39J and others. 

4.2 Folded patch antenna 

In this section a more complicated example is considered, where the farfield of a triple 
slot patch antenna fixed on a dielectric substrate is computed. The structure is taken 
from the examples coming with CST Microwave Studio as part of the CST Studio 
Suite |60|. It is illustrated in Fig. |9] with the defining points 1-13 of the patches given 
in Tab. |2] The relative permittivity of the substrate is 2.2. The substrate and metalliza- 
tion thicknesses are 0.813 mm and 0.2 mm, respectively. The example was modified 
regarding the excitation, where the original waveguide port was replaced with two dis- 
crete voltage ports, which impose a voltage across the gaps of the antenna feed at the 
position of points 2 and 3. The excitation voltage follows a Gaussian time profile with 
a standard deviation of 0.12 ns. The total simulation time is 2.5 ns. 

The farfield computation involves the determination of equivalent surface current 
densities on a collection surface F, and the subsequent solution of the Stratton-Chu 
integral under the farfield assumption, i.e., for large observation distances. We followed 
| [61j for the actual implementation and compute the equivalent current densities J(y) = 
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Table 2: Location of points 1-13 of Fig.|9]in the x-y-plane numbered from left to right 
and bottom up. 



Point 


X 1 mm 


y 1 mm 


Point 


X 1 mm 


y 1 mm 


1 


-58.5 





8 


37 


32 


2 


-1.7 





9 


-37 


34 


3 


1.7 





10 


37 


36 


4 


58.5 





11 


-37 


38 


5 


39 


30 


12 


-39 


40 


6 


-1 


32 


13 


-58.5 


60 


7 


1 


32 









n X H(y) and M(y) = E(y) x n with y a point on the collection surface and n the 
inward facing unit normal. The electric farfield is obtained by solving 



E=o(x) = ^ ^[i X M(y) + X (i X J(y))]e*'=*-^ dA, 



(40) 



where x is an observation direction, y the integration variable and k the wave num- 
ber. As time-domain simulations are performed a Fourier transform of the equivalent 
currents involving the target frequency has be carried out prior to solving ( |40] l. 

For this example, the collection surface is a box enclosing the structure at a distance 
of 2 mm. Usually the mesh is constructed such that the collection surface is obtained 
as the union of faces of a number of connected elements. In this case the elements, 
which have to be considered for solving the farfield integral can be determined in a 
preprocessing step. As this advantage cannot be exploited on adaptive meshes, we 
allow for placing the collection surface V independently of the mesh. The farfield 
integral is computed by dissecting T into (mesh independent) patches and performing 
a Gauss-Legendre quadrature on each patch. The computational domain is terminated 
by Silver-Miiller radiation boundary conditions. 



Figure 10 shows snapshots taken at a time of 0.9 ns of the electric field magnitude 
using a logarithmic color scale in the top left panel, the ft,p-mesh (top right), the el- 
ementwise error estimate (bottom left) and the element markers (bottom right). The 
viewplane is located at the bottom of the substrate. All plots show the left hand half of 
the viewplane, where the plots in the right hand panels were mirrored. In this case the 
/i-refinement level was limited to one for illustration purposes. 

Tab. [3] summarizes key figures such as the number of DoF, runtime and error esti- 
mates after the final time step for various settings. The results given in the first row (#1) 
obtained on a static mesh of third order elements are taken as the reference. The results 
#2-#4 were obtained on various adaptive meshes. Out of these, #2 is obtained using 
the same topological mesh as before and pure p-refinement with orders in between one 
to three. Results #3 and #4 were obtained on a coarser base mesh with one level of 
/i-refinement and orders of one to three and four, respectively. The error tolerance was 
set to the error of the reference solution, i.e. 0.89 • 10^^. This tolerance was met for 
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Table 3: Performance of simulations of example 4.2 using fixed and adaptive meshes 



# 


P 


L 


Elements 


DoF 

/103 


norm. 
Time 


-error 

/10-2 


TOL 
/10-2 


1 


3 


— 


11968 


4596 


1 


0.89 
(100 %) 





2 


1-3 





11968 


574- 
850 


0.15 


1.06 

(119 %) 


0.89 


3 


1-3 


0-1 


3542- 
7446 


574- 
1303 


0.27 


1.11 

(125 %) 


0.89 


4 


1-4 


0-1 


3542- 
6454 


574- 
1807 


0.34 


0.86 

(97 % ) 


0.89 


5 


2 




11968 


1939 


0.21 


1.56 

(175 %) 




6 


1 




11968 


574 


0.03 


4.60 
(517 % ) 





#4 only. Note that in the cases #2 and #3 a slightly larger global error is expected 
as the local resolution is less or at most equal compared to the reference simulation. 
The farfield is, nevertheless, very close to the reference as shown in Fig.[TT| where the 
results of #1 and #2 are shown. For adaptive simulation #4 the farfield agrees even 
better with the reference. In this context, the approach presented in | |62| is of inter- 
est, where the farfield error instead of the global solution error is employed for driving 
mesh adaptation. 



5 Conclusion 

A scheme for performing time-domain simulations with the DG method on dynamic 
/ip-adaptive meshes was proposed. The adaptation can be anisotropic in both mesh 
parameters, i.e., the element sizes h and orders p. It allows for constructing meshes, 
which meet a given error tolerance with a minimal number of degrees of freedom. The 
scheme autonomously adapts the mesh in a dynamic manner such that the error toler- 
ance is respected throughout the time-domain simulation. Apart from the reduction in 
memory and time requirements, this has the advantage that a statement about the error 
of the numerical solution can be made, therefore giving confidence in the correctness 
of the numerical result. 

The key steps of adaptive algorithms were discussed along with some popular tech- 
niques for addressing each of them. In the following an error estimate based on the idea 
of reference solutions was introduced, where the role of the reference mesh and the so- 
lution mesh was interchanged. This approach drastically reduces memory requirements 
and computing time at the cost of losing some sharpness of the estimate. Regarding el- 
ement marking a modified variable fraction strategy suitable for time-domain problems 
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was proposed. 

The attainable savings using dynamical meshes roughly scale with the multi-scale 
character of the problem at hand and can reach factors above one hundred as demon- 
strated in |27|. However, computing with adaptive meshes requires performing addi- 
tional tasks such as computing error estimates and finding the best refinement/derefinement 
option. In order to obtain a high performance scheme, the efficiency of these tasks is 
a concern. The main contribution to computational costs connected with adaptivity 
usually arise from runtime quadratures. The evaluation of the proposed error estimate 
is free of such quadratures. The same applies for projections between Finite Element 
Spaces, where precomputed integrals are employed for computing the local vectors of 
DoF during element adaptation. A number of further practical issues specific to adap- 
tive codes such as mesh and memory administration were discussed as well. Code 
profiling showed that for the presented examples the computational time consumed for 
adaptivity related tasks was around 15 % of the total computing time. Considering that 
non-adaptive simulations of comparable error levels consumed considerably more time 
(factors of three to ten in the above examples), the proposed scheme offers significant 
overall savings. 
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Figure 2: Comparison of mesh representation trees of different refinement histories. 
Starting from a single element (top) the final refinement at the bottom is obtained. For 
Option I in the left example only anisotropic refinement is applied, in Option II mixed 
anisotropic and isotropic refinement is employed. In this example, different represen- 
tation trees are obtained for identical meshes. In the right hand example identical trees 
are obtained for identical final meshes although refinements were performed in a dif- 
ferent order. Despite their identical appearance, trees III and IV differ, which becomes 
obvious when mesh derefinement is performed by cutting branches from the bottom up. 
These examples depict simple situations in two-dimensional space, in three dimensions 
more options arise. 
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Min. Tree I Min. Tree II 




Figure 3: The same refinements as in Fig. [2] are performed. The representation tree 
is constructed following a minimal depth strategy. With this strategy the maximum 
depth of the representation tree corresponds to the maximum of the refinement lev- 
els L^. Trees obtained with this strategy provide more freedom for performing mesh 
derefinement. 
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Figure 4: Illustration of the derefinement strategy using the right example of Fig. [3] 
and minimal tree III as the starting point. Only the right hand half is depicted. In 
step one the topological parent is obtained. Next, the parent is /i-refined again in order 
to generate a list of /i-candidates. In step three /ip-candidates are created by using 
different approximation orders for each /i-candidate. This step is restricted in the sense 
that all elements of one candidate have same order P. Also, each candidate has to 
have less DoF than the current mesh. All candidates are tested yielding the winning 
/ip-derefinement option. 
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Figure 5: Generation of the initial /ip-mesh for a Gauss-modulated sinusoidal wave- 
form in the fundamental mode of a rectangular waveguide using anisotropic refinement. 
The ^/-component of the electric field and the /ip-mesh is depicted in a cut view of a 
short waveguide section. The mesh is adapted iteratively such that the approximation 
error respects the tolerance TDL = 10^^. Iteration #0 shows the approximation on 
the root tesselation and the initially uniform polynomial order The adaptation termi- 
nates after 28 iterations, obtaining the initial /ip-mesh depicted at the bottom. The 
autonomous adaptation algorithm employs p-enrichment only, which is desirable as 
the solution is smooth. The /ip-meshes are depicted using a common tensor product 
visualization technique based on embedded triangles. The highest order employed 
in the initial /ip-mesh is six. The respective z-oriented edges are part of orange colored 
triangles. The maximum of Px is five (yellow). 
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Figure 6: Convergence of the global -error during construction of the initial hp-mesh 
as depicted in Fig.|5] The graph uses a logarithmic scale for the error and a linear one 
in the number of DoF, thus showing exponential convergence. The error tolerance of 
10^^ is met after 28 iterations. 
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Figure 7: Evolution of the dynamic /ip-mesh. Mesh refinement occurs predominantly 
in the form p-refinement. The snapshot show the field and mesh at the middle and end 
of the waveguide. 
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Figure 8: Temporal profiles of the global error normalized to the initial error (top, solid 
line with circles), number of elements (middle, dashed with triangles) and number of 
degrees of freedom (bottom, dotted line with diamonds). The time range covers the 
full time-domain simulation sampled at 50 instances. 




Figure 9: Triple slot folded patch antenna fixed on a dielectric substrate. The positions 
of the points in the x-y-plane are given in Tab. |2] The relative permittivity of the 
substrate is 2.2. The thicknesses of the substrate and the metallization are 0.813 mm 
and 0.2 mm, respectively. 
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Figure 10: Snapshot at time 0.9 ns of the electric field magnitude with logarithmic 
color scale (top left), the hp-mesh (top right), the estimated error (bottom left) and the 
element marker (bottom right). The viewplane is located at the bottom of the substrate. 
A non-equidistant base mesh was used for capturing the edges of the patches. Re- 
garding the adaptivity markers, an element is marked as non-refinable/irreducible only 
if no more refinement/derefinement option (h or p) is available, i.e., the given maxi- 
mum/minimum /i-level and order is met in all directions. For this demonstration the 
/i-level was limited to one and the order in between one and three. 
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Figure 11: Normalized electric farfield at a frequency of 1.5 GHz of the triple slot patch 
antenna depicted in Fig. |9] The azimuth (x-y) and elevation (z-x) plane is shown in 
the top and bottom panel. Red curves correspond to the reference solution computed 
on a non-adaptive fine mesh using third order elements (cf. Tab.[3]#l), blue curves were 
obtained with the adaptive scheme and settings according to Tab.|3]#2. 
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